suppressMessages(library(Seurat))
suppressMessages(library(ggplot2))
colorClusters <- c("#6fccc1", "#cf1fa7", "#37ee7c", "#da0322", "#ffd94a", "#207ede", "#fea111", "#b09de9",
                   "#b86fbb", "#150de2", "#ab3f2f", "#d49b6c", "#773aab", "#759474", "#e4b2b8", "#f88a40",
                   "#74b5cf", "#e6208e", "#0a3535", "#5b37c3", "#cc4025", "#d2dc7c", "#344d76", "#ba8b07",
                   "#2e10b6", "#4d0cb0", "#e1668a", "#47a58b", "#d7734a", "#ff58c7", "#edceef", "#a21d59",
                   "#8c8da5", "#ff3e5e", "#688d4b", "#214fca", "#48a8f5", "#752758", "#3b4a4a", "#11674b",
                   "#9a1bd9", "#16146d", "#7d7277", "#2d8cce", "#2460b8", "#0f9a8f", "#b11e90", "#54ed4f",
                   "#987217", "#980119", "#cf0606", "#dd182c", "#88a1e4", "#3b89ba", "#12defd", "#feba3e",
                   "#e31fa5", "#7b3537", "#098a95", "#1f1716", "#df6c94", "#9965ee", "#b438f9", "#0140f5",
                   "#e6c4a8", "#d94f7c", "#0dae56", "#b86d1d", "#0577ad", "#464551", "#6b0959", "#ccf705",
                   "#12271d")
srt <- readRDS("/mnt/nmorais-nfs/marta/pB_joana/pC_data/srt-all-samples-after-qc.rds")
srt
## An object of class Seurat 
## 16864 features across 40956 samples within 1 assay 
## Active assay: RNA (16864 features, 3000 variable features)
##  3 dimensional reductions calculated: pca, umap, tsne
table(srt$sample)
## 
## old0 old1 old3 old5  yg0  yg1  yg3  yg5 
## 6199 5596 4650 3995 5524 4434 5550 5008
DimPlot(srt, reduction = "tsne", group.by = "RNA_snn_res.0.4", cols = colorClusters) + geom_hline(yintercept= 12, color = "black", size=0.5) +
geom_hline(yintercept= 8, color = "black", size=0.5) +
geom_vline(xintercept= 7, color = "black", size=0.5) +
  geom_vline(xintercept= 15, color = "black", size=0.5) +
  ggtitle("Proliferating - group 1")

DimPlot(srt, reduction = "tsne", group.by = "RNA_snn_res.0.4", cols = colorClusters) + geom_hline(yintercept= 5, color = "black", size=0.5) +
geom_hline(yintercept= -5, color = "black", size=0.5) +
geom_vline(xintercept= 14, color = "black", size=0.5) +
  geom_vline(xintercept= 22, color = "black", size=0.5) +
  ggtitle("Proliferating - group 2")

DimPlot(srt, reduction = "tsne", group.by = "RNA_snn_res.0.4", cols = colorClusters) + geom_hline(yintercept= -8, color = "black", size=0.5) +
geom_hline(yintercept= -18, color = "black", size=0.5) +
geom_vline(xintercept= 18, color = "black", size=0.5) +
  geom_vline(xintercept= 23, color = "black", size=0.5) +
  ggtitle("Proliferating - group 3")

srt_prol <- srt[, srt$RNA_snn_res.0.4 == 12]
DimPlot(srt_prol, reduction = "tsne", group.by = "RNA_snn_res.0.4", cols = "#773aab") + geom_hline(yintercept= 12, color = "black", size=0.5) +
geom_hline(yintercept= 7, color = "black", size=0.5) +
geom_vline(xintercept= 7, color = "black", size=0.5) +
  geom_vline(xintercept= 15, color = "black", size=0.5) +
  ggtitle("Proliferating - group 1")

DimPlot(srt_prol, reduction = "tsne", group.by = "RNA_snn_res.0.4", cols = "#773aab") + geom_hline(yintercept= 5, color = "black", size=0.5) +
geom_hline(yintercept= -5, color = "black", size=0.5) +
geom_vline(xintercept= 14, color = "black", size=0.5) +
  geom_vline(xintercept= 22, color = "black", size=0.5) +
  ggtitle("Proliferating - group 2")

DimPlot(srt_prol, reduction = "tsne", group.by = "RNA_snn_res.0.4", cols = "#773aab") + geom_hline(yintercept= -8, color = "black", size=0.5) +
geom_hline(yintercept= -18, color = "black", size=0.5) +
geom_vline(xintercept= 18, color = "black", size=0.5) +
  geom_vline(xintercept= 23, color = "black", size=0.5) +
  ggtitle("Proliferating - group 3")

xcoordinates <- Embeddings(srt_prol, reduction = "tsne")[,1]
ycoordinates <- Embeddings(srt_prol, reduction = "tsne")[,2]
is.group.1 <- (xcoordinates > 7 & xcoordinates < 15) & (ycoordinates > 7 & ycoordinates < 12)
is.group.2 <- (xcoordinates > 14 & xcoordinates < 22) & (ycoordinates > -5 & ycoordinates < 5)
is.group.3 <- (xcoordinates > 18 & xcoordinates < 23) & (ycoordinates > -18 & ycoordinates < -8)
group <- as.character(srt_prol$RNA_snn_res.0.4)
group[is.group.1] <- "group 1"
group[is.group.2] <- "group 2"
group[is.group.3] <- "group 3"
group <- gsub(12, "other", group)

Number of cells per group:

table(group)
## group
## group 1 group 2 group 3   other 
##     223     403     200      28
srt_prol$group <- group
DimPlot(srt_prol, reduction = "tsne", group.by = "group") 

srt_prol <- srt_prol[, srt_prol$group != "other"]
Idents(srt_prol) <- srt_prol$group
markers <- FindAllMarkers(srt_prol,
                       test.use = "MAST")
## Calculating cluster group 3
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster group 2
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster group 1
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
markers <- markers[order(markers$avg_log2FC, decreasing = TRUE),]
markers
#saveRDS(markers, "/mnt/nmorais-nfs/marta/pB_joana/pC_data/tables-after-norm/dea-specific-groups-coordinates/proliferating-subgroup-markers.rds")
#write.table(markers, "/mnt/nmorais-nfs/marta/pB_joana/pC_data/tables-after-norm/dea-specific-groups-coordinates/proliferating-subgroup-markers.txt", append = FALSE, quote = FALSE)

Top violin plots

srt_prol$group <- factor(srt_prol$group, levels = c("group 1", "group 2", "group 3"))
Idents(srt_prol) <- srt_prol$group
VlnPlot(srt_prol, features = markers$gene[markers$cluster == "group 1"][1:10])

VlnPlot(srt_prol, features = markers$gene[markers$cluster == "group 2"][1:10])

VlnPlot(srt_prol, features = markers$gene[markers$cluster == "group 3"][1:10])

sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_3.3.6      SeuratObject_4.1.3 Seurat_4.1.0      
## 
## loaded via a namespace (and not attached):
##   [1] plyr_1.8.7                  igraph_1.3.5               
##   [3] lazyeval_0.2.2              sp_1.5-0                   
##   [5] splines_4.1.2               listenv_0.8.0              
##   [7] scattermore_0.8             GenomeInfoDb_1.30.1        
##   [9] digest_0.6.30               htmltools_0.5.3            
##  [11] fansi_1.0.3                 magrittr_2.0.3             
##  [13] tensor_1.5                  cluster_2.1.4              
##  [15] ROCR_1.0-11                 globals_0.16.1             
##  [17] matrixStats_0.62.0          spatstat.sparse_3.0-0      
##  [19] prettyunits_1.1.1           colorspace_2.0-3           
##  [21] ggrepel_0.9.1               xfun_0.34                  
##  [23] dplyr_1.0.10                crayon_1.5.2               
##  [25] RCurl_1.98-1.9              jsonlite_1.8.3             
##  [27] progressr_0.11.0            spatstat.data_3.0-0        
##  [29] survival_3.4-0              zoo_1.8-11                 
##  [31] glue_1.6.2                  polyclip_1.10-4            
##  [33] gtable_0.3.1                zlibbioc_1.40.0            
##  [35] XVector_0.34.0              leiden_0.4.3               
##  [37] DelayedArray_0.20.0         future.apply_1.9.1         
##  [39] SingleCellExperiment_1.16.0 BiocGenerics_0.40.0        
##  [41] abind_1.4-5                 scales_1.2.1               
##  [43] DBI_1.1.3                   spatstat.random_3.0-1      
##  [45] miniUI_0.1.1.1              Rcpp_1.0.9                 
##  [47] progress_1.2.2              viridisLite_0.4.1          
##  [49] xtable_1.8-4                reticulate_1.26            
##  [51] spatstat.core_2.4-4         stats4_4.1.2               
##  [53] htmlwidgets_1.5.4           httr_1.4.4                 
##  [55] RColorBrewer_1.1-3          ellipsis_0.3.2             
##  [57] ica_1.0-3                   pkgconfig_2.0.3            
##  [59] farver_2.1.1                sass_0.4.2                 
##  [61] uwot_0.1.14                 deldir_1.0-6               
##  [63] utf8_1.2.2                  tidyselect_1.2.0           
##  [65] labeling_0.4.2              rlang_1.0.6                
##  [67] reshape2_1.4.4              later_1.3.0                
##  [69] munsell_0.5.0               tools_4.1.2                
##  [71] cachem_1.0.6                cli_3.4.1                  
##  [73] generics_0.1.3              ggridges_0.5.4             
##  [75] evaluate_0.17               stringr_1.4.1              
##  [77] fastmap_1.1.0               yaml_2.3.6                 
##  [79] goftest_1.2-3               knitr_1.40                 
##  [81] fitdistrplus_1.1-8          purrr_0.3.5                
##  [83] RANN_2.6.1                  pbapply_1.5-0              
##  [85] future_1.28.0               nlme_3.1-160               
##  [87] mime_0.12                   compiler_4.1.2             
##  [89] rstudioapi_0.14             plotly_4.10.0              
##  [91] png_0.1-7                   spatstat.utils_3.0-1       
##  [93] tibble_3.1.8                bslib_0.4.0                
##  [95] stringi_1.7.8               highr_0.9                  
##  [97] lattice_0.20-45             Matrix_1.5-1               
##  [99] vctrs_0.5.0                 pillar_1.8.1               
## [101] lifecycle_1.0.3             spatstat.geom_3.0-3        
## [103] lmtest_0.9-40               jquerylib_0.1.4            
## [105] RcppAnnoy_0.0.19            data.table_1.14.4          
## [107] cowplot_1.1.1               bitops_1.0-7               
## [109] irlba_2.3.5.1               GenomicRanges_1.46.1       
## [111] httpuv_1.6.6                patchwork_1.1.2            
## [113] R6_2.5.1                    promises_1.2.0.1           
## [115] KernSmooth_2.23-20          gridExtra_2.3              
## [117] IRanges_2.28.0              parallelly_1.32.1          
## [119] codetools_0.2-18            MASS_7.3-58.1              
## [121] assertthat_0.2.1            SummarizedExperiment_1.24.0
## [123] MAST_1.20.0                 withr_2.5.0                
## [125] sctransform_0.3.5           S4Vectors_0.32.4           
## [127] GenomeInfoDbData_1.2.7      hms_1.1.2                  
## [129] mgcv_1.8-41                 parallel_4.1.2             
## [131] grid_4.1.2                  rpart_4.1.19               
## [133] tidyr_1.2.1                 rmarkdown_2.17             
## [135] MatrixGenerics_1.6.0        Rtsne_0.16                 
## [137] Biobase_2.54.0              shiny_1.7.2